%% Step 5

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Produced for JdG,AL,IM by Christopher Evans at UPF  (08/01/2018)
%
% Purpose: Allows the user to easily analyse the counterfactual exercises 
%          produced before.
%
% Output: Intended output is the level and hat variables from our algorithm
%         which can then be used in Step 6 to run. 

% NEW: no longer calls up a CounterOutput function, since not really
% needed. Also relies on Step4 baseline and results data for all variables
% + global so only call up once at begnning.
%
% Last Updated: 07/09/2019 by JDG/IM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global alphaT rhoT lambdaT etaT shockT rho sigma sigmaT lambda eta ...
    tol maxit D_hat_n Xi_hat_mnjf Xi_hat_mnj a_hat_f a_hat ...
    varphi sL_n sPi_n sD_n N J FI_J france oneminus_alpha_vector ...
    labour_share_PWT countrysecD firmsecD_sorted start start_sorted ...
    sum_by_country_dummy ISO alpha_WIOD_j alpha_WIOD_francej sector_algorithm 

%% Load data



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Homogeneity function
%
% We need to homogeneize (1-alpha)gamma given the structure of the model,
% along with the labor and intermediate shares in the data, for French
% firms. Having done this, we can then run the model as for the
% heterogeneous case.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% X_nkif: Multiplies the expenditure matrix for france with the firm dummy
%         to create a [F x N matrix]

X_nkif = firmsecD_sorted*X_mnj(start(france):start(france+1)-1,:);

% pi_X_nkif: Multiplying firm shares of expenditures (pi_nkj) by the
%            expenditure matrix that has been scaled up by the firm 
%            sector dummy

pi_X_nkif=pi_mnjf.*X_nkif;

% pi_X_nif: This is the sum over country k, should result in a [Fx1] vector

pi_X_nif=sum(pi_X_nkif,2);

%one_minus_pi_l_f: Labor shares for each firm = 1-pi_l_f

one_minus_pi_l_f=1-pi_l_f;

% Multiplying the previous summed equation by this labor share. This is
% done before multiplying by pi_M_f, which will come later.

pi_X_laborshare_nif=pi_X_nif.*one_minus_pi_l_f;

clear pi_X_nkif  X_nkif

pi_X_laborshare_nif_repmat = repmat(pi_X_laborshare_nif, [1 J*N]);     

gamma_pi_X_laborshare_nif = pi_X_laborshare_nif_repmat.*pi_M_f;

% Sum over firms. Should give a J*N matrix, which is numerator of
% homogeneous (1-alpha)*gamma

numerator = firmsecD_sorted'*gamma_pi_X_laborshare_nif;

% Sum exports over destination country per sector for france for numerator
% 32*1 matrix, which must be expanded per destination market

X_fra = X_mnj(start(france):start(france+1)-1,:);
den = sum(X_fra,2);   
denominator = repmat(den,[1 J*N]);

clear pi_X_laborshare_nif_repmat gamma_pi_X_laborshare_nif X_fra 

% Create homogeneous (1-alpha)gamma

one_minus_alp_gam_hom = numerator./denominator;

clear numerator denominator;

% Create homogeneous alpha
pi_X_laborshare_nif=pi_X_nif.*pi_l_f;
numerator = firmsecD_sorted'*pi_X_laborshare_nif;
pi_l_hom = numerator./den;

clear numerator den pi_X_laborshare_nif;
% Create homogeneous gamma given that alpha_hom = alpha_baseline at
% country-sector level

pi_l_hom_check = pi_l(start(france):start(france+1)-1,:);
corr(pi_l_hom,pi_l_hom_check) 
pi_M_hom = one_minus_alp_gam_hom./(1-repmat(pi_l_hom,[1 N*J]));

% Now fill in firm levels given homogeneous values and sectors firms in

pi_l_f_hom = firmsecD_sorted*pi_l_hom;
pi_M_f_hom = firmsecD_sorted*pi_M_hom;

clear denominator numerator one_minus_alp_gam_hom pi_l_hom pi_M_hom;

% Construct homogeneous market shares
N_jf = sum(firmsecD_sorted,1);
N_jf_repmat = repmat(N_jf,FI_J,1);
one_over_N_repmat = firmsecD_sorted ./ N_jf_repmat;
one_over_N = sum(one_over_N_repmat,2);
pi_mnjf_hom = repmat(one_over_N,1,N);

%% Rename variables
pi_mnjf0 = pi_mnjf_hom;
pi_l0 = pi_l;
pi_l_f0 = pi_l_f_hom;
pi_M_f0 = pi_M_f_hom;
pi_M0 = pi_M;
pi_c_mnj0 = pi_c_mnj;

%% MACRO VARIABLES

%X_mnj_0: Total exports from country m sector j to country n, at time t, we
%         signal this with a 0 suffix

X_mnj_0=X_mnj;

%X_mn_0 : Total exports from country m to country n (any sector). This sums
%       the sectors so we have country to country. This is at time t. 

for i=1:N
   X_mn_0(i,:)=sum(X_mnj_0(start(i):start(i+1)-1,:),1);
end

%Creating the table to show how wages in each country have changed

wage_hat=w_hat';


%X_mnj_1: This is the counterfactual expenditure matrix after
%                      adding shocks to the model. Therefore expenditure at
%                      time t+1 (after shock)
X_mnj_1=X_hat_mnj.*X_mnj_0;

%Creating the new exports from country m to country n from the
%counterfactual

for i=1:N
   X_mn_1(i,:)=nansum(X_mnj_1(start(i):start(i+1)-1,:),1);
end

X_hat_mn=X_mn_1./X_mn_0;

%X_n_0: Total exports of country m (total production), at time t

X_n_0=nansum(X_mn_0,2);

%X_mn_1: Total exports of country m after the shock, therefore time t+1

X_n_1=nansum(X_mn_1,2);

%X_hat_n: Change in total exports of country m 

X_hat_n=X_n_1./X_n_0;


% Calculating VA at the macro level:
rho_f=firmsecD_sorted*rho;
rho_nj_repmat=repmat(rho,[N 1]);
rho_mnjj_repmat=repmat(rho_nj_repmat',[N*J 1]);
rho_mnj_repmat=repmat(rho,[N N]);

pi_l_mnj_0 = repmat(pi_l0,[1 N]);
pi_l_mnj_1 = repmat(pi_l1,[1 N]);

%%
% Calculating VA (GDP)

VA_mnj_0=((1+pi_l_mnj_0.*(rho_mnj_repmat-1))./rho_mnj_repmat).*X_mnj_0;
VA_mnj_1=((1+pi_l_mnj_1.*(rho_mnj_repmat-1))./rho_mnj_repmat).*X_mnj_1;

% Summing VA over countries, to get GDP per sector 

VA_mj_0 = nansum(VA_mnj_0,2);
VA_mj_1 = nansum(VA_mnj_1,2);


% Summing the VA over sectors

for i=1:N
   VA_mn_0(i,:)=nansum(VA_mnj_0(start(i):start(i+1)-1,:),1);
   VA_mn_1(i,:)=nansum(VA_mnj_1(start(i):start(i+1)-1,:),1);
end


% VA (GDP) at the country level 

VA_n_0=nansum(VA_mn_0,2);
VA_n_1=nansum(VA_mn_1,2);

%% MICRO VARIABLES

% Firm level output 
%   - pi_mnj_1  - share of firm f in total exports from country m to n in
%   sector j. We have data when m=france. 
%  
                 
                     
% X_mnjf_0: Multiplies the expenditure matrix for france with the firm 
%           dummy to create a [F x N matrix] at time 0 (t)

X_mnjf_0 = firmsecD_sorted*X_mnj_0(start(france):start(france+1)-1,:);

% X_mnjf_1: Multiplies the expenditure matrix for france with the firm 
%           dummy to create a [F x N matrix] at time 1 (t+1)

X_mnjf_1 = firmsecD_sorted*X_mnj_1(start(france):start(france+1)-1,:);

% Multiply by firms' market shares to get sales at the firm level

X_mnjf_0=full(pi_mnjf0.*X_mnjf_0);
X_mnjf_1=full(pi_mnjf1.*X_mnjf_1);

X_hat_mnjf = X_mnjf_1./X_mnjf_0;
X_hat_mnjf(isnan(X_hat_mnjf))=1;

% 2 b) Firm f?s exports to Germany [Fx1]
% we can change this to automatically reassign the name or just select the
% correct country by specifying it in Step_5.m for now we just leave it as
% germany. Germany is the 10th country mn the WIOD

X_mDEUjf_0 = X_mnjf_0(:,10);
X_mDEUjf_1 = X_mnjf_1(:,10);

X_hat_mDEUjf = X_mDEUjf_1./X_mDEUjf_0;
X_hat_mDEUjf(isnan(X_hat_mDEUjf))=1;

% 2 c) Firm f?s ?exports? to France (i.e., domestic sales) [Fx1]

X_mFRAjf_0 = X_mnjf_0(:,france);
X_mFRAjf_1 = X_mnjf_1(:,france);

X_hat_mFRAjf = X_mFRAjf_1./X_mFRAjf_0;
X_hat_mFRAjf(isnan(X_hat_mFRAjf))=1;

% 2 d) Firm f?s exports from all other countries (aggregated) [Fx1]

X_mjf_0 = nansum(X_mnjf_0,2);
X_mjf_1 = nansum(X_mnjf_1,2);
X_hat_mjf = X_mjf_1./X_mjf_0;
X_hat_mjf(isnan(X_hat_mjf))=1;

% Value added calculation:

% pi_l_f_repmat: Replicating pi_l_f such that it spans over the N
%                        countries in our sample. 
                            
pi_l_f_repmat0=repmat(pi_l_f0, [1 N]); 
pi_l_f_repmat1=repmat(pi_l_f1, [1 N]); 
rho_fn=repmat(rho_f,[1 N]);

%Calculate VA at the firm level 

VA_fnj_0 =((1+pi_l_f_repmat0.*(rho_fn-1))./rho_fn).*X_mnjf_0;
VA_fnj_1=((1+pi_l_f_repmat1.*(rho_fn-1))./rho_fn).*X_mnjf_1; 
VA_hat_fnj=VA_fnj_1./VA_fnj_0;
VA_hat_fnj(isnan(VA_hat_fnj))=1;

% Sum across firms within a sector 
for j=1:J
    
   VA_jn_0(j,:) = nansum(VA_fnj_0(start_sorted(j):start_sorted(j+1)-1,:),1);
   VA_jn_1(j,:) = nansum(VA_fnj_1(start_sorted(j):start_sorted(j+1)-1,:),1);
   
end

% Assign the summed firm VA to country-sector level
VA_mnj_0(start(france):start(france+1)-1,:)=VA_jn_0;
VA_mnj_1(start(france):start(france+1)-1,:)=VA_jn_1;
VA_hat_mnj=VA_mnj_1./VA_mnj_0;
VA_hat_mnj(isnan(VA_hat_mnj)) = 1;

% Sum VA across destinations within a firm
VA_fj_0 = nansum(VA_fnj_0,2);
VA_fj_1 = nansum(VA_fnj_1,2);
VA_hat_fj=VA_fj_1./VA_fj_0;

% Sum across firms within a sector
for j=1:J
    VA_mj_0(start(france)+j-1,:)=nansum(VA_fj_0(start_sorted(j):start_sorted(j+1)-1,:),1);
    VA_mj_1(start(france)+j-1,:)=nansum(VA_fj_1(start_sorted(j):start_sorted(j+1)-1,:),1); 
end
VA_hat_mj= VA_mj_1./VA_mj_0;
VA_hat_mj(isnan(VA_hat_mj)) = 1;

% Assign french VA to country-country level 
VA_francen_0 = nansum(VA_fnj_0,1);
VA_francen_1 = nansum(VA_fnj_1,1);
VA_mn_0(france,:)=VA_francen_0;
VA_mn_1(france,:)=VA_francen_1;   

% Assign french VA to country level 
VA_francen_0 = nansum(VA_fj_0,1);
VA_francen_1 = nansum(VA_fj_1,1);
VA_n_0(france)=VA_francen_0;
VA_n_1(france)=VA_francen_1;
VA_hat_n=VA_n_1./VA_n_0;

          
%% Labor market compensation:


% calculating labour market compensation:


% real_wage_income_hat: This is taken from notes 2 section 1.3, the
% equation on total labor compensation.

% Calculating wage income at the country*sector level from the baseline:
wage_income_mj_0 = ...
            (pi_l0.*(rho_nj_repmat-1)./rho_nj_repmat).*sum(X_mnj_0,2);

% Aggregating across sectors
for i=1:N
    wage_income_n_0(i) = sum(wage_income_mj_0(start(i):start(i+1)-1),1);
end

% Calculating wage income at the country level from the counterfactual:

wage_income_mj_1 = ...
            (pi_l1.*(rho_nj_repmat-1)./rho_nj_repmat).*sum(X_mnj_1,2);   

        
for i=1:N
 wage_income_n_1(i) = sum(wage_income_mj_1(start(i):start(i+1)-1),1);
end        

wage_income_hat_n = wage_income_n_1./wage_income_n_0;

% Real wage income hat for country m is total wage income hat over price
% hat

real_wage_income_cpi_hat_n =(wage_income_hat_n./P_hat_n)';

%% Double GDP Deflator calculations
[ DoubleDeflation_deflator_n,RGDP_doubledef_hat_n, RGDP_doubledef_hat_fj] = ...
    DoubleDeflation_deflatorCES_fun_approx(P_hat_mnj,...
        pi_M0,pi_M_f0,pi_l_f0,pi_l_f1,pi_l0,pi_l1,...
        wage_hat,X_mnj_1,X_mnj_0,X_mjf_0,X_mjf_1,VA_n_0,VA_hat_fj,VA_hat_n,rho_f,rho_nj_repmat);
real_wage_doubledef_hat = wage_hat./DoubleDeflation_deflator_n;
real_wage_income_doubledef_hat_n =(wage_income_hat_n'./DoubleDeflation_deflator_n);

%% Real GDP Deflator calculations

[ GDP_deflator_n,RGDP_def_hat_n, RGDP_def_hat_fj] = GDP_deflatorCES_fun_approx(P_hat_mnj,...
    pi_M_f0,pi_l_f0,wage_hat,pi_M0,pi_l0,X_hat_mnj,...
    VA_n_0,VA_mnj_0,VA_hat_n,X_hat_mnjf,VA_fnj_0,VA_hat_fj);
real_wage_def_hat = wage_hat./GDP_deflator_n;
real_wage_income_def_hat_n =(wage_income_hat_n'./GDP_deflator_n);

%% Real variables CPI Deflator

real_wage_cpi_hat = wage_hat./P_hat_n';
RGDP_cpi_hat_n=VA_hat_n./P_hat_n';
RGDP_cpi_hat_fj =VA_hat_fj./P_hat_n(france);


% Calculating Price hat and Real GDP by country-sector pair

sigma_mj_repmat=repmat(sigma,[N 1]);
sigma_mnj_repmat=repmat(sigma_mj_repmat,[1 N]);
  
% Calculating P_hat_nj 

% P_hat_pi_c: Multiply the price hat with the share of final good
%             consumption 

P_hat_pi_c= P_hat_mnj.^(1-sigma_mnj_repmat).*pi_c_mnj0;

% sum_P_hat_pi_c: Summing the above equation by country m 

sum_P_hat_pi_c=sum_by_country_dummy*P_hat_pi_c;

% P_hat_nj: Taking the expression to the power of 1/(1-rho). It gives us
%           the new Price hat but for country n and sector j

P_hat_nj= sum_P_hat_pi_c.^(1./(1-repmat(sigma,[1 N])));

% Set the value =1 if it is equal to infinity. The infinity problem comes 
% as we are putting 0^(negative number), the zero is due to no final goods
% purchased in that sector and country. 

P_hat_nj(P_hat_nj==inf)=1;

P_hat_nj_france=P_hat_nj(:,france);

for i=1:N
    VA_hat_mj_matrix(:,i) = VA_hat_mj(start(i):start(i+1)-1,1);
end
RGDP_cpi_hat_mj = VA_hat_mj_matrix ./P_hat_nj;
RGDP_cpi_hat_mj_france = RGDP_cpi_hat_mj(:,france);


%% Openess Calculations
%turning off domestic channel

X_mnj_dom_0 = X_mnj_0;

for i=1:N
    X_mnj_dom_0(start(i):start(i+1)-1,i)=0;
end

% Bilateral Exports
for i=1:N
   X_mn_dom_0(i,:)=nansum(X_mnj_dom_0(start(i):start(i+1)-1,:),1);
end

% Aggregate exports
X_n_dom_0=nansum(X_mn_dom_0,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Exports at country m sector j (summed across all foreign
% destinations)

X_mj_dom_0=nansum(X_mnj_dom_0,2);

% Imports at country n sector j at post-wedge matrix 
% at the sector level 

M_nj_dom_0= sum_by_country_dummy*X_mnj_dom_0;

% Imports by country n , for the sector model 

M_n_dom_0=nansum(M_nj_dom_0,1);

% Exports by all for the sector model 
exports_over_GDP_all=X_n_dom_0./VA_n_0;
exports_over_GDP_france=exports_over_GDP_all(france);

% Imports by all for the sector model 
imports_over_GDP_all=M_n_dom_0'./VA_n_0;
imports_over_GDP_france=imports_over_GDP_all(france);

% Calculating exports over GDP for france at the sector level using
% sectoral exports and GDP by sector 

exports_over_GDP_france_sector=...
    X_mj_dom_0(start(france):start(france+1)-1)./...
    VA_mj_0(start(france):start(france+1)-1);


% Imports over GDP for france at the secotr levle using sectoral imports
% and GDP by sector

imports_over_GDP_france_sector= ...
    M_nj_dom_0(:,france)./...
    VA_mj_0(start(france):start(france+1)-1);

% Sanity test:

sum_VA_f_0=nansum(VA_fj_0);

sum_VA_f_1=nansum(VA_fj_1);

sum_VA_hat_f=sum_VA_f_1./sum_VA_f_0 ;

%This should line up with real french gdp (it does)
sum_VA_hat_f_real_cpi=sum_VA_hat_f./P_hat_n(france) 
RGDP_cpi_hat_n(france)

clear sum_VA_hat_f sum_VA_f_1 sum_VA_f_0;

%% TFP Calculation

% Change in labour supply is given by: 
% L hat = (wage hat/ price hat)^(1/(psi-1) 

L_hat_n = (real_wage_cpi_hat).^(1/(varphi-1));

% TFP_hat_n: Change in TFP, assuming Y = TFP * L ; 
TFP_hat_cpi_n= RGDP_cpi_hat_n./L_hat_n;

% TFP_hat_france: isolating the TFP for France
TFP_hat_cpi_france=TFP_hat_cpi_n(france);

% Calculating the TFP using the GDP deflator as the price hat instead.

% Change in labour supply is given by: 
% L hat = (wage hat/ price hat)^(1/(psi-1) same for deflator or cpi

L_hat_deflator_n=L_hat_n;

% TFP_hat_deflator_n: Calculating TFP as Y_hat/L_hat = TFP_hat
TFP_hat_def_n = RGDP_def_hat_n./L_hat_deflator_n;

% TFP_hat_deflator_france: Isolating TFP for france
TFP_hat_def_france = TFP_hat_def_n(france);

% Calculating the TFP using the GDP double-deflator as the price hat instead.

% Change in labour supply is given by: 
% L hat = (wage hat/ price hat)^(1/(psi-1) same for deflator or cpi

L_hat_doubledef_n=L_hat_n;

% TFP_hat_deflator_n: Calculating TFP as Y_hat/L_hat = TFP_hat
TFP_hat_doubledef_n = RGDP_doubledef_hat_n./L_hat_doubledef_n;

% TFP_hat_deflator_france: Isolating TFP for france
TFP_hat_doubledef_france = TFP_hat_doubledef_n(france);

% Calculating Labour share at country-sector level, using september 24th
% draft, bottom of page 12, the labour composition:

% pi_l_mnjf_0: labour share multiplied by expenditure X at the
% firm level

pi_l_mnjf_0=pi_l_f_repmat0.*X_mnjf_0;
pi_l_mnjf_1=pi_l_f_repmat1.*X_mnjf_1;

% sum_france_pi_l_X_mnj_0: Summing over the firms 

sum_france_pi_l_X_mnj_0 = firmsecD_sorted'*pi_l_mnjf_0;
sum_france_pi_l_X_mnj_1 = firmsecD_sorted'*pi_l_mnjf_1;

% sum_pi_l_X_mnj_0: Summing over other counries k

sum_pi_l_X_mnj_0= nansum(sum_france_pi_l_X_mnj_0,2);
sum_pi_l_X_mnj_1= nansum(sum_france_pi_l_X_mnj_1,2);


% wage_L_hat_france_j: The hat (change pre and post shock) is wage
% multiplied by labour at the country-sector level

wage_L_hat_france_j =sum_pi_l_X_mnj_1./sum_pi_l_X_mnj_0;


% L_hat_france_j: Labour supply to each sector, found dividing the above
% expression by the wage in France

L_hat_france_j = wage_L_hat_france_j./wage_hat(france);


%TFP_france_j: TFP in france by sector is Y_hat_sector/L_hat_sector

TFP_hat_cpi_france_j = RGDP_cpi_hat_mj_france./L_hat_france_j;
